Stan for Asymptotic Recovery Curves

We can use non-linear mixed effects models for normally distributed responses to estimate asymptotic recovery curves for VIQ and PIQ separately. It’s feasible to extend these models to handle VIQ and PIQ as a multivariate response, but doing so is not straighforward.

Here, we use Stan to build a model for a multivariate response which will allow us to compare the parameters for the two models in a way that exploits the covariance between the responses.

We start with univariate models which we then combine into a multivariate model.

Loading packages:

library(spida2)
library(magrittr, pos = 1000) 
library(lattice)
library(latticeExtra)

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

Explore data

data(iq)
?iq
head(iq)
  DAYSPC DCOMA  SEX    AgeIn   ID PIQ VIQ
1     30     4 Male 20.67077 3358  87  89
2     16    17 Male 55.28816 3535  95  77
3     40     1 Male 55.91513 3547  95 116
4     13    10 Male 61.66461 3592  59  73
5     19     6 Male 30.12731 3728  67  73
6     13     3 Male 57.06229 3790  76  69
names(iq) <- tolower(names(iq))

dd <-iq

(p <- xyplot(piq ~ dayspc | sex, dd, groups = id, type = 'b'))

update(p, xlim = c(0,4000))

if(interactive) {
  rows <- numeric(0)
  # can repeat:
  trellis.focus()
  rows <- c(rows, panel.identify(labels=dd$id))
  # end
  trellis.unfocus()
  rows
  save(rows, file = 'rows.rda')
}

load('rows.rda', verbose = T)
Loading objects:
  rows
dd[rows,] %>% sortdf(~dcoma+dayspc)
    dayspc dcoma  sex    agein   id piq viq
47    3333     9 Male 43.93977 2600  86  80
327   3337     9 Male 43.93977 2600 101  84
303   1491    21 Male 22.00684  651  71  94
325   3412    21 Male 22.00684  651  68  92
310   1926   130 Male 28.27379 1939  95 108
321   3111   130 Male 28.27379 1939  88 111
326   3864   130 Male 28.27379 1939  88 105
# id = 2600 retested 4 days apart
dd[dd$id %in% dd[rows,'id'],]  %>%  sortdf(~dcoma+dayspc)
    dayspc dcoma  sex    agein   id piq viq
47    3333     9 Male 43.93977 2600  86  80
327   3337     9 Male 43.93977 2600 101  84
303   1491    21 Male 22.00684  651  71  94
325   3412    21 Male 22.00684  651  68  92
234    295   130 Male 28.27379 1939  67 117
271    562   130 Male 28.27379 1939  85 111
310   1926   130 Male 28.27379 1939  95 108
321   3111   130 Male 28.27379 1939  88 111
326   3864   130 Male 28.27379 1939  88 105
if(interactive){
  library(p3d)
  Init3d()
  dd$dcoma.cat <- cut(dd$dcoma, c(-1,2,5,10,20,50,Inf))
  Plot3d( viq ~ piq + log(dayspc) |
            dcoma.cat, dd, groups = id,
          col = heat.colors(6))
  Plot3d( viq ~ log(dcoma+2) + log(dayspc) |
            dcoma.cat, dd, groups = id,
          col = heat.colors(6))
  Plot3d( piq ~ log(dcoma+2) + log(dayspc) |
            dcoma.cat, dd, groups = id,
          col = heat.colors(12)[1:6])
  fg()
  Id3d()
}

Stan program

stanfile <- 'asymp_model_1.stan' 
pr <- function(x) print(readLines(x), quote = FALSE) # write a function for repetitive tasks
pr(stanfile)
 [1]                                                         
 [2] // Asymptotic recovery model                            
 [3] // - random intercept                                   
 [4]                                                         
 [5] data {                                                  
 [6]   int N;                                                
 [7]   int J;                                                
 [8]   vector[N] y;                                          
 [9]   vector[N] time;                                       
[10]   vector[N] coma;                                       
[11]   int id[N];                                            
[12]   int Np;  // predictor values for prediction           
[13]   vector[Np] time_p;                                    
[14]   vector[Np] coma_p;                                    
[15] }                                                       
[16] transformed data{                                       
[17]   real ln2;                                             
[18]   ln2 = log(2); // factor for half-recovery time        
[19] }                                                       
[20] parameters {                                            
[21]   real hrt;    // half-recovery time                    
[22]   real asymp;  // mean asymptotic recovery if dcoma = 0 
[23]   real bcoma;  // duration of coma parameter            
[24]   real init_def; // initial deficit                     
[25]   vector[J] u;   // random intercept residual           
[26]   real <lower=0> sigma;   // sd within                  
[27]   real <lower=0> sigma_u; // sd between                 
[28] }                                                       
[29] transformed parameters {                                
[30]   real estd_reliability;                                
[31]   estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[32] }                                                       
[33] model {                                                 
[34]   u ~ normal(0,sigma_u);                                
[35]   y ~ normal(asymp + u[id] + bcoma * coma +             
[36]              init_def * exp(-time/(hrt*ln2)), sigma);   
[37] }                                                       
[38] generated quantities {                                  
[39]   vector[Np] y_fit;                                     
[40]   y_fit = asymp + bcoma * coma_p +                      
[41]              init_def * exp(-time_p/(hrt*ln2));         
[42] }                                                       
  system.time(
  asymp_model_dso <- stan_model(stanfile)
)
   user  system elapsed 
  0.094   0.035   0.842 

Prepare data list

names(dd)
[1] "dayspc" "dcoma"  "sex"    "agein"  "id"     "piq"    "viq"   
pred <- expand.grid(time=seq(0, 3*365, 30), 
                    coma = seq(0,30,5))
dat <- list(
  N = nrow(dd),
  id = nid <- as.numeric(as.factor(dd$id)),
  J = max(nid),
  y = dd$piq,
  time = dd$dayspc,
  coma = sqrt(dd$dcoma),
  Np = nrow(pred),
  time_p = pred$time,
  coma_p = sqrt(pred$coma)
)
set.seed(789723)
mod <- sampling(asymp_model_dso, dat, chains = 4, iter = 2000)
# Note: chains and iter are defaults
# traceplot(mod)
names(mod)
  [1] "hrt"              "asymp"            "bcoma"           
  [4] "init_def"         "u[1]"             "u[2]"            
  [7] "u[3]"             "u[4]"             "u[5]"            
 [10] "u[6]"             "u[7]"             "u[8]"            
 [13] "u[9]"             "u[10]"            "u[11]"           
 [16] "u[12]"            "u[13]"            "u[14]"           
 [19] "u[15]"            "u[16]"            "u[17]"           
 [22] "u[18]"            "u[19]"            "u[20]"           
 [25] "u[21]"            "u[22]"            "u[23]"           
 [28] "u[24]"            "u[25]"            "u[26]"           
 [31] "u[27]"            "u[28]"            "u[29]"           
 [34] "u[30]"            "u[31]"            "u[32]"           
 [37] "u[33]"            "u[34]"            "u[35]"           
 [40] "u[36]"            "u[37]"            "u[38]"           
 [43] "u[39]"            "u[40]"            "u[41]"           
 [46] "u[42]"            "u[43]"            "u[44]"           
 [49] "u[45]"            "u[46]"            "u[47]"           
 [52] "u[48]"            "u[49]"            "u[50]"           
 [55] "u[51]"            "u[52]"            "u[53]"           
 [58] "u[54]"            "u[55]"            "u[56]"           
 [61] "u[57]"            "u[58]"            "u[59]"           
 [64] "u[60]"            "u[61]"            "u[62]"           
 [67] "u[63]"            "u[64]"            "u[65]"           
 [70] "u[66]"            "u[67]"            "u[68]"           
 [73] "u[69]"            "u[70]"            "u[71]"           
 [76] "u[72]"            "u[73]"            "u[74]"           
 [79] "u[75]"            "u[76]"            "u[77]"           
 [82] "u[78]"            "u[79]"            "u[80]"           
 [85] "u[81]"            "u[82]"            "u[83]"           
 [88] "u[84]"            "u[85]"            "u[86]"           
 [91] "u[87]"            "u[88]"            "u[89]"           
 [94] "u[90]"            "u[91]"            "u[92]"           
 [97] "u[93]"            "u[94]"            "u[95]"           
[100] "u[96]"            "u[97]"            "u[98]"           
[103] "u[99]"            "u[100]"           "u[101]"          
[106] "u[102]"           "u[103]"           "u[104]"          
[109] "u[105]"           "u[106]"           "u[107]"          
[112] "u[108]"           "u[109]"           "u[110]"          
[115] "u[111]"           "u[112]"           "u[113]"          
[118] "u[114]"           "u[115]"           "u[116]"          
[121] "u[117]"           "u[118]"           "u[119]"          
[124] "u[120]"           "u[121]"           "u[122]"          
[127] "u[123]"           "u[124]"           "u[125]"          
[130] "u[126]"           "u[127]"           "u[128]"          
[133] "u[129]"           "u[130]"           "u[131]"          
[136] "u[132]"           "u[133]"           "u[134]"          
[139] "u[135]"           "u[136]"           "u[137]"          
[142] "u[138]"           "u[139]"           "u[140]"          
[145] "u[141]"           "u[142]"           "u[143]"          
[148] "u[144]"           "u[145]"           "u[146]"          
[151] "u[147]"           "u[148]"           "u[149]"          
[154] "u[150]"           "u[151]"           "u[152]"          
[157] "u[153]"           "u[154]"           "u[155]"          
[160] "u[156]"           "u[157]"           "u[158]"          
[163] "u[159]"           "u[160]"           "u[161]"          
[166] "u[162]"           "u[163]"           "u[164]"          
[169] "u[165]"           "u[166]"           "u[167]"          
[172] "u[168]"           "u[169]"           "u[170]"          
[175] "u[171]"           "u[172]"           "u[173]"          
[178] "u[174]"           "u[175]"           "u[176]"          
[181] "u[177]"           "u[178]"           "u[179]"          
[184] "u[180]"           "u[181]"           "u[182]"          
[187] "u[183]"           "u[184]"           "u[185]"          
[190] "u[186]"           "u[187]"           "u[188]"          
[193] "u[189]"           "u[190]"           "u[191]"          
[196] "u[192]"           "u[193]"           "u[194]"          
[199] "u[195]"           "u[196]"           "u[197]"          
[202] "u[198]"           "u[199]"           "u[200]"          
[205] "sigma"            "sigma_u"          "estd_reliability"
[208] "y_fit[1]"         "y_fit[2]"         "y_fit[3]"        
[211] "y_fit[4]"         "y_fit[5]"         "y_fit[6]"        
[214] "y_fit[7]"         "y_fit[8]"         "y_fit[9]"        
[217] "y_fit[10]"        "y_fit[11]"        "y_fit[12]"       
[220] "y_fit[13]"        "y_fit[14]"        "y_fit[15]"       
[223] "y_fit[16]"        "y_fit[17]"        "y_fit[18]"       
[226] "y_fit[19]"        "y_fit[20]"        "y_fit[21]"       
[229] "y_fit[22]"        "y_fit[23]"        "y_fit[24]"       
[232] "y_fit[25]"        "y_fit[26]"        "y_fit[27]"       
[235] "y_fit[28]"        "y_fit[29]"        "y_fit[30]"       
[238] "y_fit[31]"        "y_fit[32]"        "y_fit[33]"       
[241] "y_fit[34]"        "y_fit[35]"        "y_fit[36]"       
[244] "y_fit[37]"        "y_fit[38]"        "y_fit[39]"       
[247] "y_fit[40]"        "y_fit[41]"        "y_fit[42]"       
[250] "y_fit[43]"        "y_fit[44]"        "y_fit[45]"       
[253] "y_fit[46]"        "y_fit[47]"        "y_fit[48]"       
[256] "y_fit[49]"        "y_fit[50]"        "y_fit[51]"       
[259] "y_fit[52]"        "y_fit[53]"        "y_fit[54]"       
[262] "y_fit[55]"        "y_fit[56]"        "y_fit[57]"       
[265] "y_fit[58]"        "y_fit[59]"        "y_fit[60]"       
[268] "y_fit[61]"        "y_fit[62]"        "y_fit[63]"       
[271] "y_fit[64]"        "y_fit[65]"        "y_fit[66]"       
[274] "y_fit[67]"        "y_fit[68]"        "y_fit[69]"       
[277] "y_fit[70]"        "y_fit[71]"        "y_fit[72]"       
[280] "y_fit[73]"        "y_fit[74]"        "y_fit[75]"       
[283] "y_fit[76]"        "y_fit[77]"        "y_fit[78]"       
[286] "y_fit[79]"        "y_fit[80]"        "y_fit[81]"       
[289] "y_fit[82]"        "y_fit[83]"        "y_fit[84]"       
[292] "y_fit[85]"        "y_fit[86]"        "y_fit[87]"       
[295] "y_fit[88]"        "y_fit[89]"        "y_fit[90]"       
[298] "y_fit[91]"        "y_fit[92]"        "y_fit[93]"       
[301] "y_fit[94]"        "y_fit[95]"        "y_fit[96]"       
[304] "y_fit[97]"        "y_fit[98]"        "y_fit[99]"       
[307] "y_fit[100]"       "y_fit[101]"       "y_fit[102]"      
[310] "y_fit[103]"       "y_fit[104]"       "y_fit[105]"      
[313] "y_fit[106]"       "y_fit[107]"       "y_fit[108]"      
[316] "y_fit[109]"       "y_fit[110]"       "y_fit[111]"      
[319] "y_fit[112]"       "y_fit[113]"       "y_fit[114]"      
[322] "y_fit[115]"       "y_fit[116]"       "y_fit[117]"      
[325] "y_fit[118]"       "y_fit[119]"       "y_fit[120]"      
[328] "y_fit[121]"       "y_fit[122]"       "y_fit[123]"      
[331] "y_fit[124]"       "y_fit[125]"       "y_fit[126]"      
[334] "y_fit[127]"       "y_fit[128]"       "y_fit[129]"      
[337] "y_fit[130]"       "y_fit[131]"       "y_fit[132]"      
[340] "y_fit[133]"       "y_fit[134]"       "y_fit[135]"      
[343] "y_fit[136]"       "y_fit[137]"       "y_fit[138]"      
[346] "y_fit[139]"       "y_fit[140]"       "y_fit[141]"      
[349] "y_fit[142]"       "y_fit[143]"       "y_fit[144]"      
[352] "y_fit[145]"       "y_fit[146]"       "y_fit[147]"      
[355] "y_fit[148]"       "y_fit[149]"       "y_fit[150]"      
[358] "y_fit[151]"       "y_fit[152]"       "y_fit[153]"      
[361] "y_fit[154]"       "y_fit[155]"       "y_fit[156]"      
[364] "y_fit[157]"       "y_fit[158]"       "y_fit[159]"      
[367] "y_fit[160]"       "y_fit[161]"       "y_fit[162]"      
[370] "y_fit[163]"       "y_fit[164]"       "y_fit[165]"      
[373] "y_fit[166]"       "y_fit[167]"       "y_fit[168]"      
[376] "y_fit[169]"       "y_fit[170]"       "y_fit[171]"      
[379] "y_fit[172]"       "y_fit[173]"       "y_fit[174]"      
[382] "y_fit[175]"       "y_fit[176]"       "y_fit[177]"      
[385] "y_fit[178]"       "y_fit[179]"       "y_fit[180]"      
[388] "y_fit[181]"       "y_fit[182]"       "y_fit[183]"      
[391] "y_fit[184]"       "y_fit[185]"       "y_fit[186]"      
[394] "y_fit[187]"       "y_fit[188]"       "y_fit[189]"      
[397] "y_fit[190]"       "y_fit[191]"       "y_fit[192]"      
[400] "y_fit[193]"       "y_fit[194]"       "y_fit[195]"      
[403] "y_fit[196]"       "y_fit[197]"       "y_fit[198]"      
[406] "y_fit[199]"       "y_fit[200]"       "y_fit[201]"      
[409] "y_fit[202]"       "y_fit[203]"       "y_fit[204]"      
[412] "y_fit[205]"       "y_fit[206]"       "y_fit[207]"      
[415] "y_fit[208]"       "y_fit[209]"       "y_fit[210]"      
[418] "y_fit[211]"       "y_fit[212]"       "y_fit[213]"      
[421] "y_fit[214]"       "y_fit[215]"       "y_fit[216]"      
[424] "y_fit[217]"       "y_fit[218]"       "y_fit[219]"      
[427] "y_fit[220]"       "y_fit[221]"       "y_fit[222]"      
[430] "y_fit[223]"       "y_fit[224]"       "y_fit[225]"      
[433] "y_fit[226]"       "y_fit[227]"       "y_fit[228]"      
[436] "y_fit[229]"       "y_fit[230]"       "y_fit[231]"      
[439] "y_fit[232]"       "y_fit[233]"       "y_fit[234]"      
[442] "y_fit[235]"       "y_fit[236]"       "y_fit[237]"      
[445] "y_fit[238]"       "y_fit[239]"       "y_fit[240]"      
[448] "y_fit[241]"       "y_fit[242]"       "y_fit[243]"      
[451] "y_fit[244]"       "y_fit[245]"       "y_fit[246]"      
[454] "y_fit[247]"       "y_fit[248]"       "y_fit[249]"      
[457] "y_fit[250]"       "y_fit[251]"       "y_fit[252]"      
[460] "y_fit[253]"       "y_fit[254]"       "y_fit[255]"      
[463] "y_fit[256]"       "y_fit[257]"       "y_fit[258]"      
[466] "y_fit[259]"       "lp__"            
pars <- grepv('^u|y_fit', names(mod), invert = T)
pars
[1] "hrt"              "asymp"            "bcoma"           
[4] "init_def"         "sigma"            "sigma_u"         
[7] "estd_reliability" "lp__"            
traceplot(mod, pars = pars)

pairs(mod, pars = pars)

if(interactive) {
  library(shinystan)
  mod_sso <- launch_shinystan(mod)
}
zd <- as.data.frame(mod, pars = pars)
head(zd)
       hrt     asymp     bcoma  init_def    sigma  sigma_u
1 352.0651 102.14958 -2.409008 -17.19386 7.255366 11.40271
2 280.1005 100.73336 -2.460842 -17.74048 7.549275 14.95228
3 202.0059  96.23615 -1.543932 -18.50024 6.701158 12.54657
4 254.4917 100.87878 -1.832296 -20.22619 6.789999 12.59399
5 284.2113 100.09366 -2.042092 -17.07062 7.634592 13.40865
6 247.9092  96.78811 -1.456427 -17.42886 6.868731 12.95669
  estd_reliability      lp__
1        0.7118160 -1426.493
2        0.7968666 -1443.540
3        0.7780491 -1416.420
4        0.7747862 -1414.997
5        0.7551779 -1417.888
6        0.7806170 -1401.731
xyplot(hrt~asymp, zd)

print(mod, pars = pars)
Inference for Stan model: asymp_model_1.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                     mean se_mean    sd     2.5%      25%      50%
hrt                277.27    1.07 64.93   173.72   231.85   267.65
asymp              100.96    0.05  1.99    97.13    99.63   100.92
bcoma               -1.92    0.01  0.40    -2.69    -2.20    -1.93
init_def           -19.50    0.04  1.87   -23.25   -20.73   -19.50
sigma                6.92    0.01  0.44     6.13     6.61     6.90
sigma_u             13.19    0.01  0.84    11.59    12.62    13.16
estd_reliability     0.78    0.00  0.03     0.71     0.76     0.78
lp__             -1415.25    0.52 14.77 -1444.76 -1425.33 -1414.83
                      75%    97.5% n_eff Rhat
hrt                315.79   430.61  3692 1.00
asymp              102.29   104.94  1567 1.00
bcoma               -1.64    -1.13  1342 1.00
init_def           -18.27   -15.75  2800 1.00
sigma                7.21     7.84  1486 1.00
sigma_u             13.74    14.95  4000 1.00
estd_reliability     0.81     0.84  1834 1.00
lp__             -1404.74 -1388.23   810 1.01

Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:47:50 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
get_posterior_mean(mod, pars = pars)[,5]
             hrt            asymp            bcoma         init_def 
      277.267594       100.955620        -1.920296       -19.503431 
           sigma          sigma_u estd_reliability             lp__ 
        6.921886        13.193070         0.782237     -1415.249592 

Note that you can do extensive plotting of the posterior sample with ‘shinystan’ including some not highly satisfying 3d plots.

Random initial deficit

It wasn’t feasible to fit a model that included random initial deficits using ‘lme’ so this will be an interesting experiment.

Note that the data will be identical.

stanfile <- 'asymp_model_2.stan' 
pr(stanfile)
 [1]                                                         
 [2] // Asymptotic recovery model                            
 [3] // - random intercept                                   
 [4]                                                         
 [5] data {                                                  
 [6]   int N;                                                
 [7]   int J;                                                
 [8]   vector[N] y;                                          
 [9]   vector[N] time;                                       
[10]   vector[N] coma;                                       
[11]   int id[N];                                            
[12]   int Np;  // predictor values for prediction           
[13]   vector[Np] time_p;                                    
[14]   vector[Np] coma_p;                                    
[15] }                                                       
[16] transformed data{                                       
[17]   real ln2;                                             
[18]   ln2 = log(2); // factor for half-recovery time        
[19] }                                                       
[20] parameters {                                            
[21]   real hrt;    // half-recovery time                    
[22]   real asymp;  // mean asymptotic recovery if dcoma = 0 
[23]   real bcoma;  // duration of coma parameter            
[24]   real init_def_mean; // initial deficit                
[25]   vector[J] u;   // random intercept residual           
[26]   vector[J] u_init_def; // random initial deficit       
[27]   real <lower=0> sigma;   // sd within                  
[28]   real <lower=0> sigma_u; // sd between                 
[29]   real <lower=0> sigma_idef;                            
[30] }                                                       
[31] transformed parameters {                                
[32]   real estd_reliability;                                
[33]   estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[34] }                                                       
[35] model {                                                 
[36]   // note few changes from non-random init_def          
[37]   u ~ normal(0,sigma_u);                                
[38]   // add an empirical prior for u_init_def:             
[39]   u_init_def ~ normal(init_def_mean, sigma_idef);       
[40]   // add u_init_def[id] and element-wise multiplication 
[41]   // of vectors (.*):                                   
[42]   y ~ normal(asymp + u[id] + bcoma * coma +             
[43]              (init_def_mean + u_init_def[id])           
[44]                 .* exp(-time/(hrt*ln2)), sigma);        
[45] }                                                       
[46] generated quantities {                                  
[47]   vector[Np] y_fit;                                     
[48]   y_fit = asymp + bcoma * coma_p +                      
[49]              init_def_mean * exp(-time_p/(hrt*ln2));    
[50] }                                                       
system.time(
  asymp_model_2_dso <- stan_model(stanfile)
)
In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                 from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file385317a558f2.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition
   user  system elapsed 
 33.396   2.446  36.522 

Fitting the model. The data is the same.

mod2 <- sampling(asymp_model_2_dso, dat, chains = 4, iter = 2000)
Warning: There were 745 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
pars <- grepv('^u|y_fit', names(mod2), invert = T)
pars
[1] "hrt"              "asymp"            "bcoma"           
[4] "init_def_mean"    "sigma"            "sigma_u"         
[7] "sigma_idef"       "estd_reliability" "lp__"            
traceplot(mod2, pars = pars)

pairs(mod2, pars = pars)

Poor convergence: possible remedies

The variability in the sd of initial deficit is poorly estimated.

We can try to improve this in a least two ways:

  1. Reparametrize the model: Taking the initial point at time 0 puts it outside the range of the data. Using a later time as time 0 could improve the estimate. This can be done just by changing the data and subtracting a constant from time. We then need to remember the changed interpretation of parameters although most parameters are not affected, e.g. ‘asymp’ and ‘hrt’
  2. Use an informative prior for ‘sigma_idef’. So far all models have used flat priors for hyper-parameters. Gelman would not hesitate to use much more informative priors.

We’ll try step 1 because it’s easy. We’ll subtract 180 from time so that time 0 will be relocated to approximately 6 months.

Reparametrize

dat_6 <- list(
  N = nrow(dd),
  id = nid <- as.numeric(as.factor(dd$id)),
  J = max(nid),
  y = dd$piq,
  time = dd$dayspc - 180,
  coma = sqrt(dd$dcoma),
  Np = nrow(pred),
  time_p = pred$time - 180,
  coma_p = sqrt(pred$coma)
)

mod2_6 <- sampling(asymp_model_2_dso, dat_6, chains = 4, iter = 2000)
Warning: There were 815 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
pars <- grepv('^u|y_fit', names(mod2_6), invert = T)
pars
[1] "hrt"              "asymp"            "bcoma"           
[4] "init_def_mean"    "sigma"            "sigma_u"         
[7] "sigma_idef"       "estd_reliability" "lp__"            
traceplot(mod2_6, pars = pars)

pairs(mod2_6, pars = pars)

Informative prior

We will try a Gamma prior on ‘sigma_idef’

stanfile <- 'asymp_model_3.stan' 
pr(stanfile)
 [1]                                                         
 [2] // Asymptotic recovery model                            
 [3] // - random intercept                                   
 [4]                                                         
 [5] data {                                                  
 [6]   int N;                                                
 [7]   int J;                                                
 [8]   vector[N] y;                                          
 [9]   vector[N] time;                                       
[10]   vector[N] coma;                                       
[11]   int id[N];                                            
[12]   int Np;  // predictor values for prediction           
[13]   vector[Np] time_p;                                    
[14]   vector[Np] coma_p;                                    
[15]   // gamma prior parameters                             
[16]   real gamma_df;                                        
[17]   real gamma_scale;                                     
[18] }                                                       
[19] transformed data{                                       
[20]   real ln2;                                             
[21]   ln2 = log(2); // factor for half-recovery time        
[22] }                                                       
[23] parameters {                                            
[24]   real hrt;    // half-recovery time                    
[25]   real asymp;  // mean asymptotic recovery if dcoma = 0 
[26]   real bcoma;  // duration of coma parameter            
[27]   real init_def_mean; // initial deficit                
[28]   vector[J] u;   // random intercept residual           
[29]   vector[J] u_init_def; // random initial deficit       
[30]   real <lower=0> sigma;   // sd within                  
[31]   real <lower=0> sigma_u; // sd between                 
[32]   real <lower=0> sigma_idef;                            
[33] }                                                       
[34] transformed parameters {                                
[35]   real estd_reliability;                                
[36]   estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;
[37] }                                                       
[38] model {                                                 
[39]   // Gamma prior                                        
[40]   sigma_idef ~ gamma(gamma_df, 1/gamma_scale);          
[41]   // note few changes from non-random init_def          
[42]   u ~ normal(0,sigma_u);                                
[43]   // add an empirical prior for u_init_def:             
[44]   u_init_def ~ normal(init_def_mean, sigma_idef);       
[45]   // add u_init_def[id] and element-wise multiplication 
[46]   // of vectors (.*):                                   
[47]   y ~ normal(asymp + u[id] + bcoma * coma +             
[48]              (init_def_mean + u_init_def[id])           
[49]                 .* exp(-time/(hrt*ln2)), sigma);        
[50] }                                                       
[51] generated quantities {                                  
[52]   vector[Np] y_fit;                                     
[53]   y_fit = asymp + bcoma * coma_p +                      
[54]              init_def_mean * exp(-time_p/(hrt*ln2));    
[55] }                                                       
system.time(
  asymp_model_3_dso <- stan_model(stanfile)
)
In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                 from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file385352cf9003.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition
   user  system elapsed 
 33.660   0.914  35.301 

Fitting the model. The data is the same plus the two parameters for the gamma prior.

mod3_6 <- sampling(asymp_model_3_dso, 
                 c(dat_6, gamma_df = 4, gamma_scale = .5), 
                 chains = 4, iter = 2000,
                 control = list(adapt_delta = 0.85))
Warning: There were 87 divergent transitions after warmup. Increasing adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Warning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
pars <- grepv('^u|y_fit', names(mod3_6), invert = T)
pars
[1] "hrt"              "asymp"            "bcoma"           
[4] "init_def_mean"    "sigma"            "sigma_u"         
[7] "sigma_idef"       "estd_reliability" "lp__"            
traceplot(mod3_6, pars = pars)

pairs(mod3_6, pars = pars)
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

print(mod3_6, pars = pars)
Inference for Stan model: asymp_model_3.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                         mean      se_mean           sd          2.5%
hrt               1.56740e+02 1.032100e+02 1.521100e+02  6.750000e+00
asymp             5.09500e+01 3.539000e+01 5.008000e+01  1.700000e-01
bcoma            -9.20000e-01 6.900000e-01 1.010000e+00 -2.540000e+00
init_def_mean    -2.62000e+00 9.700000e-01 1.460000e+00 -5.300000e+00
sigma             3.75000e+00 2.220000e+00 3.160000e+00  1.400000e-01
sigma_u           8.87000e+00 3.210000e+00 4.570000e+00  2.410000e+00
sigma_idef        1.60000e+00 8.700000e-01 1.280000e+00  2.500000e-01
estd_reliability  8.80000e-01 7.000000e-02 1.000000e-01  7.200000e-01
lp__             -1.38287e+32 1.554271e+32 4.849026e+32 -2.718517e+33
                           25%           50%      75%    97.5% n_eff  Rhat
hrt               1.448000e+01  7.911000e+01   290.33   424.99     2  3.33
asymp             1.260000e+00  4.822000e+01   100.99   104.18     2 38.71
bcoma            -1.890000e+00 -2.700000e-01    -0.01     0.14     2  3.68
init_def_mean    -3.910000e+00 -1.640000e+00    -1.43    -1.00     2  2.81
sigma             8.700000e-01  3.330000e+00     6.85     7.65     2 10.26
sigma_u           5.720000e+00  8.720000e+00    13.11    14.52     2  8.29
sigma_idef        2.500000e-01  1.240000e+00     3.58     3.58     2  3.73
estd_reliability  7.900000e-01  9.300000e-01     0.98     1.00     2  4.42
lp__             -2.679679e+31 -3.924177e+11 -1556.21 -1362.73    10  1.38

Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:54:15 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
get_posterior_mean(mod3_6, pars = pars)[,5] %>% cbind
                             .
hrt               1.567368e+02
asymp             5.094795e+01
bcoma            -9.159970e-01
init_def_mean    -2.620795e+00
sigma             3.747419e+00
sigma_u           8.870717e+00
sigma_idef        1.603489e+00
estd_reliability  8.841998e-01
lp__             -1.382870e+32
get_posterior_mean(mod2_6, pars = pars)[,5] %>% cbind
                             .
hrt               1.778507e+04
asymp             9.013245e+01
bcoma            -1.981526e+00
init_def_mean    -1.050906e+01
sigma             6.558797e+00
sigma_u           8.059216e+00
sigma_idef        2.496417e+00
estd_reliability  4.656921e-01
lp__             -3.148898e+31
get_posterior_mean(mod2, pars = pars)[,5] %>% cbind
                             .
hrt                276.8901869
asymp              100.6649413
bcoma               -1.8387664
init_def_mean       -9.7724863
sigma                6.8567078
sigma_u             13.1091332
sigma_idef           3.7960164
estd_reliability     0.7827252
lp__             -1733.0739153

Baysian uniform prior

We’ll just use our intuition about IQs. The standard deviation couldn’t ‘possibly’ be less than 3 and couldn’t ‘possibly’ be greater than 30!

stanfile <- 'asymp_model_3b.stan' 
pr(stanfile)
 [1]                                                                                    
 [2] // Asymptotic recovery model                                                       
 [3] // - random intercept                                                              
 [4] // - 'Bayesian' prior on sigma_idef                                                
 [5]                                                                                    
 [6] data {                                                                             
 [7]   int N;                                                                           
 [8]   int J;                                                                           
 [9]   vector[N] y;                                                                     
[10]   vector[N] time;                                                                  
[11]   vector[N] coma;                                                                  
[12]   int id[N];                                                                       
[13]   int Np;  // predictor values for prediction                                      
[14]   vector[Np] time_p;                                                               
[15]   vector[Np] coma_p;                                                               
[16] }                                                                                  
[17] transformed data{                                                                  
[18]   real ln2;                                                                        
[19]   ln2 = log(2); // factor for half-recovery time                                   
[20] }                                                                                  
[21] parameters {                                                                       
[22]   real hrt;    // half-recovery time                                               
[23]   real asymp;  // mean asymptotic recovery if dcoma = 0                            
[24]   real bcoma;  // duration of coma parameter                                       
[25]   real init_def_mean; // initial deficit                                           
[26]   vector[J] u;   // random intercept residual                                      
[27]   vector[J] u_init_def; // random initial deficit                                  
[28]   real <lower=0> sigma;   // sd within                                             
[29]   real <lower=0> sigma_u; // sd between                                            
[30]   real <lower=3,upper=30> sigma_idef; // uniform prior from 3 to 30                
[31] }                                                                                  
[32] transformed parameters {                                                           
[33]   real estd_reliability;                                                           
[34]   estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;                           
[35] }                                                                                  
[36] model {                                                                            
[37]   // Gamma prior                                                                   
[38]   // sigma_idef ~ gamma(gamma_df, 1/gamma_scale); // commenting out the gamma prior
[39]   // note few changes from non-random init_def                                     
[40]   u ~ normal(0,sigma_u);                                                           
[41]   // add an empirical prior for u_init_def:                                        
[42]   u_init_def ~ normal(init_def_mean, sigma_idef);                                  
[43]   // add u_init_def[id] and element-wise multiplication                            
[44]   // of vectors (.*):                                                              
[45]   y ~ normal(asymp + u[id] + bcoma * coma +                                        
[46]              (init_def_mean + u_init_def[id])                                      
[47]                 .* exp(-time/(hrt*ln2)), sigma);                                   
[48] }                                                                                  
[49] generated quantities {                                                             
[50]   vector[Np] y_fit;                                                                
[51]   y_fit = asymp + bcoma * coma_p +                                                 
[52]              init_def_mean * exp(-time_p/(hrt*ln2));                               
[53] }                                                                                  
system.time(
  asymp_model_3b_dso <- stan_model(stanfile)
)
In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                 from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file385363a7af4e.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition
   user  system elapsed 
 33.790   0.897  35.409 

Fitting the model. The data is the same plus the two parameters for the gamma prior.

mod3b_6 <- sampling(asymp_model_3b_dso, 
                   dat_6,  
                   chains = 4, iter = 2000,
                   control = list(adapt_delta = 0.85))
Warning: There were 691 divergent transitions after warmup. Increasing adapt_delta above 0.85 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
pars <- grepv('^u|y_fit', names(mod3b_6), invert = T)
pars
[1] "hrt"              "asymp"            "bcoma"           
[4] "init_def_mean"    "sigma"            "sigma_u"         
[7] "sigma_idef"       "estd_reliability" "lp__"            
traceplot(mod3b_6, pars = pars)

pairs(mod3b_6, pars = pars)
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'
Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

Warning in KernSmooth::bkde2D(x, bandwidth = bandwidth, gridsize = nbin, :
Binning grid too coarse for current (small) bandwidth: consider increasing
'gridsize'

print(mod3b_6, pars = pars)
Inference for Stan model: asymp_model_3b.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                          mean      se_mean           sd          2.5%
hrt               3.356300e+02 2.641500e+02 6.938500e+02  2.590000e+00
asymp             5.138000e+01 3.605000e+01 5.102000e+01  1.500000e-01
bcoma            -1.280000e+00 6.700000e-01 9.800000e-01 -2.520000e+00
init_def_mean    -1.960000e+00 2.590000e+00 3.720000e+00 -7.610000e+00
sigma             1.889000e+01 1.668000e+01 2.768000e+01  4.900000e-01
sigma_u           7.170000e+00 3.550000e+00 5.070000e+00  1.400000e-01
sigma_idef        8.900000e+00 3.460000e+00 4.990000e+00  3.050000e+00
estd_reliability  6.000000e-01 2.600000e-01 3.700000e-01  0.000000e+00
lp__             -6.936327e+80 7.594519e+80 1.612219e+81 -6.538357e+81
                           25%      50%      75%    97.5% n_eff  Rhat
hrt               4.153000e+01   137.78   522.59  1038.61     7  1.22
asymp             5.200000e-01    48.31   102.34   106.14     2 37.69
bcoma            -1.920000e+00    -1.75    -0.14     0.43     2  3.83
init_def_mean    -5.520000e+00    -0.68     1.68     1.68     2  5.84
sigma             5.000000e-01     7.14     8.66    79.36     3  2.74
sigma_u           3.270000e+00     4.37    12.10    14.20     2  7.41
sigma_idef        3.900000e+00     9.40    13.29    15.81     2  7.76
estd_reliability  1.400000e-01     0.73     0.91     0.99     2  9.36
lp__             -1.156439e+80 -2078.99 -1795.19 -1710.55     5  2.32

Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:56:50 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
get_posterior_mean(mod3b_6, pars = pars)[,5] %>% cbind
                             .
hrt               3.356334e+02
asymp             5.137587e+01
bcoma            -1.284750e+00
init_def_mean    -1.959387e+00
sigma             1.889447e+01
sigma_u           7.166087e+00
sigma_idef        8.904122e+00
estd_reliability  6.036895e-01
lp__             -6.936327e+80
get_posterior_mean(mod2_6, pars = pars)[,5] %>% cbind
                             .
hrt               1.778507e+04
asymp             9.013245e+01
bcoma            -1.981526e+00
init_def_mean    -1.050906e+01
sigma             6.558797e+00
sigma_u           8.059216e+00
sigma_idef        2.496417e+00
estd_reliability  4.656921e-01
lp__             -3.148898e+31
get_posterior_mean(mod2, pars = pars)[,5] %>% cbind
                             .
hrt                276.8901869
asymp              100.6649413
bcoma               -1.8387664
init_def_mean       -9.7724863
sigma                6.8567078
sigma_u             13.1091332
sigma_idef           3.7960164
estd_reliability     0.7827252
lp__             -1733.0739153

Alas, this didn’t work either.

Sensitivity analysis

Since it looks like we can’t estimate ‘sigma_idef’ from the data with our model but we are concerned that it might have an uncertain impact on our other results, we can resort to a sensitivity analysis. Try a variety of plausible values and see how it affects other results.

We also put a plausible lower bound of 10 on sigma_u.

stanfile <- 'asymp_model_3c.stan' 
pr(stanfile)
 [1]                                                                                        
 [2] // Asymptotic recovery model                                                           
 [3] // - random intercept                                                                  
 [4] // - sensitivity analysis for sigma_idef                                               
 [5]                                                                                        
 [6] data {                                                                                 
 [7]   int N;                                                                               
 [8]   int J;                                                                               
 [9]   vector[N] y;                                                                         
[10]   vector[N] time;                                                                      
[11]   vector[N] coma;                                                                      
[12]   int <lower=1,upper=J> id[N];                                                         
[13]   int Np;  // predictor values for prediction                                          
[14]   vector[Np] time_p;                                                                   
[15]   vector[Np] coma_p;                                                                   
[16]   real sigma_idef;                                                                     
[17] }                                                                                      
[18] transformed data{                                                                      
[19]   real ln2;                                                                            
[20]   ln2 = log(2); // factor for half-recovery time                                       
[21] }                                                                                      
[22] parameters {                                                                           
[23]   real hrt;    // half-recovery time                                                   
[24]   real asymp;  // mean asymptotic recovery if dcoma = 0                                
[25]   real bcoma;  // duration of coma parameter                                           
[26]   real init_def_mean; // initial deficit                                               
[27]   vector[J] u;   // random intercept residual                                          
[28]   vector[J] u_init_def; // random initial deficit                                      
[29]   real <lower=0.0> sigma;   // sd within                                               
[30]   real <lower=10.0> sigma_u; // sd between  // lower limit on sigma_u                  
[31]   // real <lower=3,upper=30> sigma_idef; // will input as data for sensitivity analysis
[32] }                                                                                      
[33] transformed parameters {                                                               
[34]   real estd_reliability;                                                               
[35]   estd_reliability = sigma_u^2 / (sigma_u^2 + sigma^2) ;                               
[36] }                                                                                      
[37] model {                                                                                
[38]   // Gamma prior                                                                       
[39]   // sigma_idef ~ gamma(gamma_df, 1/gamma_scale); // commenting out the gamma prior    
[40]   // note few changes from non-random init_def                                         
[41]   u ~ normal(0,sigma_u);                                                               
[42]   // add an empirical prior for u_init_def:                                            
[43]   u_init_def ~ normal(init_def_mean, sigma_idef);                                      
[44]   // add u_init_def[id] and element-wise multiplication                                
[45]   // of vectors (.*):                                                                  
[46]   y ~ normal(asymp + u[id] + bcoma * coma +                                            
[47]              (init_def_mean + u_init_def[id])                                          
[48]                 .* exp(-time/(hrt*ln2)), sigma);                                       
[49] }                                                                                      
[50] generated quantities {                                                                 
[51]   vector[Np] y_fit;                                                                    
[52]   y_fit = asymp + bcoma * coma_p +                                                     
[53]              init_def_mean * exp(-time_p/(hrt*ln2));                                   
[54] }                                                                                      
system.time(
  asymp_model_3c_dso <- stan_model(stanfile)
)
   user  system elapsed 
  0.098   0.030   0.837 

Fitting the model. The data is the same plus the two parameters for the gamma prior.

mod3c_6 <- sampling(asymp_model_3c_dso, 
                    c(dat_6, sigma_idef = .1),  
                    chains = 4, iter = 2000)
Warning: There were 72 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
pars <- grepv('^u|y_fit', names(mod3c_6), invert = T)
pars
[1] "hrt"              "asymp"            "bcoma"           
[4] "init_def_mean"    "sigma"            "sigma_u"         
[7] "estd_reliability" "lp__"            
traceplot(mod3c_6, pars = pars)

pairs(mod3c_6, pars = pars)

print(mod3c_6, pars = pars)
Inference for Stan model: asymp_model_3c.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                          mean      se_mean           sd          2.5%
hrt               2.186900e+02 8.997000e+01 1.358900e+02  2.740000e+00
asymp             7.633000e+01 3.039000e+01 4.301000e+01  1.900000e+00
bcoma            -1.680000e+00 3.400000e-01 5.900000e-01 -2.730000e+00
init_def_mean    -2.560000e+00 1.610000e+00 2.340000e+00 -5.140000e+00
sigma             5.320000e+00 1.980000e+00 2.830000e+00  4.700000e-01
sigma_u           1.301000e+01 2.100000e-01 7.600000e-01  1.173000e+01
estd_reliability  8.400000e-01 7.000000e-02 1.000000e-01  7.100000e-01
lp__             -4.106383e+76 4.787314e+76 8.285904e+76 -2.918382e+77
                           25%      50%      75%    97.5% n_eff  Rhat
hrt               1.020900e+02   259.77   313.00   417.63     2  2.68
asymp             7.105000e+01   100.31   101.92   104.76     2 27.26
bcoma            -2.130000e+00    -1.75    -0.89    -0.89     3  1.67
init_def_mean    -4.220000e+00    -3.56    -0.73     1.33     2  3.89
sigma             4.380000e+00     6.72     7.12     7.81     2  7.69
sigma_u           1.254000e+01    12.81    13.51    14.73    13  1.08
estd_reliability  7.700000e-01     0.80     0.91     1.00     2  3.46
lp__             -1.159598e+76 -1524.22 -1508.49 -1484.01     3  3.38

Samples were drawn using NUTS(diag_e) at Wed Mar 29 13:59:04 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
get_posterior_mean(mod3c_6, pars = pars)[,5] %>% cbind
                             .
hrt               2.186874e+02
asymp             7.633317e+01
bcoma            -1.680734e+00
init_def_mean    -2.564241e+00
sigma             5.323000e+00
sigma_u           1.301404e+01
estd_reliability  8.352306e-01
lp__             -4.106383e+76
get_posterior_mean(mod3c_6, pars = pars)[,5] %>% cbind
                             .
hrt               2.186874e+02
asymp             7.633317e+01
bcoma            -1.680734e+00
init_def_mean    -2.564241e+00
sigma             5.323000e+00
sigma_u           1.301404e+01
estd_reliability  8.352306e-01
lp__             -4.106383e+76

Multivariate model for VIQ and PIQ

Instead of simple standard deviations for each VIQ and PIQ, we need covariance matrices to describe the within- and the between-subject covariance between VIQ and PIQ.

The first model uses a uniform prior on covariance matrices using constraints generated by Stan on covariance matrices.

It is more efficient to use a combination of an LKJ prior on correlation together with priors of variances. The second model below use this form of parametrization.

stanfile <- 'asymp_model_4.stan' 
pr(stanfile)
 [1] //                                                                     
 [2] //  Multivariate model for VIQ and PIQ                                 
 [3] //                                                                     
 [4]                                                                        
 [5] data {                                                                 
 [6]   int N;                                                               
 [7]   int J;                                                               
 [8]   matrix[N,2] iq;                                                      
 [9]   vector[N] time;                                                      
[10]   vector[N] coma;                                                      
[11]   int id[N];                                                           
[12] }                                                                      
[13] transformed data{                                                      
[14]   real ln2;                                                            
[15]   vector[2] zero;                                                      
[16]   ln2 = log(2);                                                        
[17]   for ( i in 1:2) zero[i] = 0.0;                                       
[18] }                                                                      
[19] parameters {                                                           
[20]   vector <lower=1,upper=10000>[2] hrt;                                 
[21]   vector <lower=0,upper=200>[2] asymp;                                 
[22]   vector <lower=-100,upper=100>[2] init_def;                           
[23]   vector [2] bcoma;                                                    
[24]   vector[2] u[J];                                                      
[25]   cov_matrix[2] Sigma;                                                 
[26]   cov_matrix[2] Sigma_u;                                               
[27] }                                                                      
[28] transformed parameters {                                               
[29]   real hrt_diff;                                                       
[30]   real bcoma_diff;                                                     
[31]   hrt_diff = hrt[2] - hrt[1];                                          
[32]   bcoma_diff = bcoma[2] - bcoma[1];                                    
[33] }                                                                      
[34] model {                                                                
[35]   vector[2] eta;                                                       
[36]   // for the multinormal distribution we need to loop over observations
[37]   for(j in 1:J) u[j] ~ multi_normal(zero, Sigma_u);                    
[38]   for(n in 1:N) {                                                      
[39]     eta[1] = asymp[1] + u[id[n],1] + bcoma[1] * coma[n] +              
[40]             init_def[1] * exp(-time[n]/(hrt[1]*ln2));                  
[41]     eta[2] = asymp[2] + u[id[n],2] + bcoma[2] * coma[n] +              
[42]             init_def[2] * exp(-time[n]/(hrt[2]*ln2));                  
[43]     iq[n,] ~ multi_normal(eta, Sigma);                                 
[44]   }                                                                    
[45] }                                                                      
[46]                                                                        
system.time(
  asymp_model_4_dso <- stan_model(stanfile)
)
In file included from /usr/local/lib/R/site-library/BH/include/boost/config.hpp:39:0,
                 from /usr/local/lib/R/site-library/BH/include/boost/math/tools/config.hpp:13,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math.hpp:4,
                 from /usr/local/lib/R/site-library/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file38533f792847.cpp:8:
/usr/local/lib/R/site-library/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined [enabled by default]
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition
   user  system elapsed 
 36.302   1.233  38.270 

Data for multivariate model

dat <- list(
  N = nrow(dd),
  id = nid <- as.numeric(as.factor(dd$id)),
  J = max(nid),
  iq = cbind(dd$viq,dd$piq),
  time = dd$dayspc,
  coma = sqrt(dd$dcoma)
)

mult_mod <- sampling( asymp_model_4_dso, dat)

pars <- grepv('^u' ,names(mult_mod), invert = T)
traceplot(mult_mod, pars = pars)

pairs(mult_mod, pars = pars)
the following parameters were dropped because they are duplicative
Sigma[1,2] Sigma_u[1,2]

print(mult_mod, pars = pars)
Inference for Stan model: asymp_model_4.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                 mean se_mean    sd     2.5%      25%      50%      75%
hrt[1]          66.00    0.45 19.35    36.99    52.62    63.39    76.24
hrt[2]         246.89    0.94 51.88   161.55   209.98   242.37   278.12
asymp[1]        99.69    0.06  1.60    96.70    98.58    99.65   100.75
asymp[2]       100.35    0.07  1.94    96.58    99.06   100.33   101.64
init_def[1]    -23.23    0.11  4.74   -34.41   -25.74   -22.62   -19.94
init_def[2]    -19.43    0.05  2.02   -23.14   -20.78   -19.43   -18.12
bcoma[1]        -0.68    0.02  0.39    -1.45    -0.94    -0.68    -0.43
bcoma[2]        -1.91    0.02  0.42    -2.74    -2.19    -1.91    -1.63
Sigma[1,1]      33.10    0.11  4.35    25.73    29.97    32.70    35.88
Sigma[2,1]      20.33    0.12  4.29    12.73    17.39    20.02    22.80
Sigma[1,2]      20.33    0.12  4.29    12.73    17.39    20.02    22.80
Sigma[2,2]      49.60    0.18  6.66    38.24    45.06    48.87    53.62
Sigma_u[1,1]   162.71    0.39 19.54   128.64   149.15   161.34   175.20
Sigma_u[2,1]   119.30    0.38 17.84    87.32   107.00   118.07   130.66
Sigma_u[1,2]   119.30    0.38 17.84    87.32   107.00   118.07   130.66
Sigma_u[2,2]   175.87    0.46 22.68   135.33   160.06   174.46   190.14
hrt_diff       180.89    0.80 50.81    98.17   145.43   175.23   210.50
bcoma_diff      -1.22    0.01  0.33    -1.85    -1.45    -1.23    -1.00
lp__         -2603.56    0.71 21.04 -2645.33 -2618.05 -2602.95 -2589.04
                97.5% n_eff Rhat
hrt[1]         111.44  1864 1.00
hrt[2]         362.33  3063 1.00
asymp[1]       102.93   620 1.01
asymp[2]       104.26   804 1.00
init_def[1]    -15.79  1844 1.00
init_def[2]    -15.30  1658 1.00
bcoma[1]         0.11   657 1.01
bcoma[2]        -1.08   774 1.00
Sigma[1,1]      42.72  1495 1.00
Sigma[2,1]      29.84  1256 1.00
Sigma[1,2]      29.84  1256 1.00
Sigma[2,2]      64.49  1394 1.00
Sigma_u[1,1]   203.72  2518 1.00
Sigma_u[2,1]   156.86  2226 1.00
Sigma_u[1,2]   156.86  2226 1.00
Sigma_u[2,2]   224.79  2440 1.00
hrt_diff       295.17  4000 1.00
bcoma_diff      -0.58  3180 1.00
lp__         -2564.57   877 1.00

Samples were drawn using NUTS(diag_e) at Wed Mar 29 14:05:08 2017.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Using a Cholesky parametrization